3. 1D Cartesian structured VDS along parallel lines

In this notebook, we perform pseudo-random variable density sampling along the phase encoding direction. A handcrafted density is designed and samples are then drawn by virtually inverting its cumulative density function. Then these samples define the selected phase encoding lines retained in the sampling mask. Low frequencies are more sampled than higher frequencies.

• Author: Philippe Ciuciu (philippe.ciuciu@cea.fr)
• Date: 06/24/2022
• Date : 04/02/2025 (use of BrainWeb database using multi-contrast T1/T2/etc imaging)
• Target: IEEE EMBS-SPS Summer School on Novel acquisition and image reconstruction strategies in accelerated Magnetic Resonance Imaging

#DISPLAY BRAIN PHANTOM
'''
在医学成像（如 MRI、CT）和信号处理中，幻影（Phantom）是指用于测试、校准或研究的人工数据或物理模型。
它通常是一个合成（模拟）图像或实验装置，用来模拟真实世界的结构，比如大脑、心脏或其他人体组织。
'''
# 1. 导入必要的库

%matplotlib inline

import numpy as np
import os.path as op
import os
import math ; import cmath
import matplotlib.pyplot as plt
import sys

from skimage import data, io, filters       # 用于图像处理。

import pywt as pw
import matplotlib.pyplot as plt

import brainweb_dl as bwdl                  # 用于下载和处理脑部 MRI 幻影数据。

# 2. 设置 Matplotlib 显示参数
plt.rcParams["image.origin"]="lower"        # 设置图像的原点在左下角（默认是左上角）。
plt.rcParams["image.cmap"]='Greys_r'        # 设置图像颜色映射为灰度反转（黑白）。

# 3. 下载并处理 MRI 图像
mri_img = bwdl.get_mri(4, "T1")[70, ...].astype(np.float32)
'''
bwdl.get_mri(4, "T1") : 从 BrainWeb 数据库下载一个 T1 加权 MRI 扫描
70 : 获取 MRI 体数据（3D）中的第 70 层切片，即 2D 切片图像。
.astype(np.float32) : 转换数据类型为 float32，以便进行计算。
''' 
#mri_img = bwdl.get_mri(4, "T2")[120, ...].astype(np.float32)

# 4. 显示 MRI 图像
print(mri_img.shape)                        # (256, 256)
img_size = mri_img.shape[0]                 # 获取图像的宽度（假设是正方形图像）。
plt.figure()

# 5. 绘制 MRI 图像
plt.imshow(abs(mri_img))                    # 显示 MRI 图像，abs() 主要用于防止负像素值影响显示。
plt.title("Original brain image")           
plt.show()

(256, 256)

import numpy as np
import sys
import numpy.random as ra

def Reconstruction (mask, signoise, decay, nb_samples):

    '''
    mask == "VDS_along_parallel_lines"
    signoise = 10
    decay = 1.
    nb_samples = (int)(img_size/4): 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
    '''

    # -1. 定义快速傅里叶变换 (FFT) 和 逆快速傅里叶变换 (IFFT)
    #import numpy.fft as fft
    norm = "ortho"
    #norm = None
    def fft(x):
        return np.fft.fft2(x, norm=norm)

    def ifft(x):
        return np.fft.ifft2(x, norm=norm)
    
    if mask == "VDS_along_parallel_lines" :
        eps = sys.float_info.epsilon 

        # 0. 初始化一个 k-space 采样掩码矩阵
        c = np.ones((1, img_size), dtype="float64")     
        kspace_mask = np.tile(c, (img_size, 1)) 


        # 1. 生成伪随机的采样密度和 CDF（累计分布函数）
        kspace_lines = np.linspace(-1/2., 1/2.,img_size)*img_size   
        # define the taret sampling density (ie non-iniform over k-space lines)

        ## decay = 1.

        # Define the sampling density
        p_decay = np.power(np.abs(kspace_lines),-decay)     
        p_decay = p_decay/np.sum(p_decay)
        # generate its CDF
        cdf_pdecay = np.cumsum(p_decay)

        pmax = p_decay.max()
        pmin = p_decay.min()

        '''
        # Plot the density and its cumulative distribution function (CDF) p_decay
        plt.figure()
        fig, axs = plt.subplots(1, 2, figsize= (5,5) )
        axs[0].plot(kspace_lines, p_decay)
        axs[1].plot(kspace_lines, cdf_pdecay)
        '''

        # 2. 生成伪随机采样线 
        # Perform pseudo-random sampling: technique used:
        # draw uniform variables and invert the CDF to get back to p_decay-distributed sampled
        ## nb_samples = (int)(img_size/4)            # 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
        '''
        print(nb_samples)
        '''
        samples = ra.uniform(0, 1, nb_samples)      # 在 [0, 1] 范围内生成 nb_samples 个均匀分布的随机数。 
                                                    # 每个随机数的值都在 0 到 1 之间，且每个随机数的概率是相等的。
        gen_klines = [int(kspace_lines[np.argwhere(cdf_pdecay == min(cdf_pdecay[(cdf_pdecay - r) > 0]))]) for r in samples]

        # 3. 对采样线进行偏移和平移
        # shift the samples lines by half of the k-space (ie image size as we're in Cartesian ref)
        '''
        生成的采样线 gen_klines 在 [-img_size/2, img_size/2] 范围内，
        需要平移到 [0, img_size] 以匹配图像矩阵的索引格式。
        '''
        gen_klinesb = ((np.array(gen_klines) - 1) / 1).astype(int) + (int)(img_size/2)
        #gen_klinesb = ((np.array(gen_klines) - 1) / 1).astype(int)
        '''
        print(gen_klinesb)   # unsorted samples
        '''

        '''
        # 4. 绘制采样密度分布和直方图
        times = np.arange(1, img_size, 1)
        # np.bincount(gen_klinesb) 计算每条采样线的出现次数，得到伪随机采样线的实际分布。
        lc = np.bincount(gen_klinesb, minlength=len(times))     

        # 绘制采样线分布直方图并与原始密度分布 p_decay 比较，检查采样分布是否与目标分布一致。
        # check that histogram of sample values fits the prescribed density p_decay
        plt.figure()
        plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Assigned k-space lines')
        plot2, = plt.plot(p_decay,'g',label='Original Spectrum')
        plt.xlabel('k-lines')
        plt.ylabel('Probability')
        plt.legend(handles=[plot1,plot2])
        plt.show()
        '''

        # 5. 生成采样掩码 (Sampling Mask) : kspace_mask
        #print(p_decay.min())
        sampled_klines = np.array(np.unique(gen_klinesb))       # np.unique() 用于去除重复值，返回唯一的采样线索引列表。
        '''
        print(sampled_klines)
        '''

        nblines = np.size(sampled_klines)
        '''
        print(nblines)
        '''
        threshold = 2. * p_decay.min()  # sys.float_info.epsilon \simeq 2e-16
        kspace_mask = np.zeros((img_size,img_size), dtype="float64")
        kspace_mask[sampled_klines,:] = np.ones((nblines,img_size) , dtype="float64")

        '''
        plt.figure()
        plt.imshow(kspace_mask, cmap='gray')
        plt.show()
        '''

        # 6. 添加复值高斯噪声
        # Generate the kspace data: first Fourier transform the image
        kspace_data = np.fft.fftshift(fft(mri_img))
        #add Gaussian complex-valued random noise
        ## signoise = 10
        #kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
        # Simulate independent noise realization on the real & imag parts
        kspace_data += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise

        # 7. 使用采样掩码进行次采样
        # Mask data to perform subsampling
        kspace_data *= kspace_mask

        # 8. 进行 MRI 重建
        # Zero order solution
        image_rec0 = ifft(np.fft.ifftshift(kspace_data))

        return kspace_data, image_rec0
                                                    
    

$\leadsto$   Here we write a function using a mask based on pseudo-random variable density sampling along the phase encoding direction, derived from handcrafted densities. The samples are then drawn by virtually inverting its cumulative density function (CDF). Then these samples define the selected phase encoding lines retained in the sampling mask.

The key step is to generate sample lines gen_klines within [-img_size/2, img_size/2] :

gen_klines = [int(kspace_lines[np.argwhere(cdf_pdecay == min(cdf_pdecay[(cdf_pdecay - r) > 0]))]) for r in samples] :

For each random sample r in samples, we find the first element in cdf_pdecay that is greater than r , then retrieve the corresponding line position from kspace_lines and and convert it to the integer format int .

• where cdf_pdecay is obtained by computing the probability density function (PDF) p_decay = np.power(np.abs(kspace_lines),-decay) , followed by normalization as p_decay = p_decay/np.sum(p_decay) , integrating this PDF cdf_pdecay = np.cumsum(p_decay)

• samples is a set of nb_samples uniformly distributed random numbers generated within the range [0, 1] : samples = ra.uniform(0, 1, nb_samples)

We shift the samples lines by half of the k-space (i.e., image size as we're in Cartesian ref). It indicates where to put the 1 in the sampling mask kspace_mask initialized by 0.

We apply this mask to the k-space data for subsampling. Finally, the reconstructed image is obtained by applying the inverse Fourier transform to the masked k-space data.

Q1  

• Second block: Vary parameter decay in the definition of the sampling density and study its impact on image quality (SSIM score).

mask="VDS_along_parallel_lines"
vector_signoise = [0, 10, 20, 100]
vector_decay = [0.1, 0.5, 0.8, 1, 1.5, 4., 7.]
vector_nb_samples = [int(img_size / i) for i in [2, 4, 8, 16, 32]]

'''
mask == "VDS_along_parallel_lines"
signoise = 10
decay = 1.
nb_samples = (int)(img_size/4): 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
'''
fig, axs = plt.subplots(len(vector_decay), 2, figsize=(24, 24))

for i in range(len(vector_decay)) :
    for j in range(2):

        signoise = vector_signoise[1]
        decay = vector_decay[i]
        nb_samples = vector_nb_samples[1]

        kspace_data, image_rec0 = Reconstruction (mask, signoise, decay, nb_samples)

        # 显示被噪声影响的 k-space 数据
        axs[i, 0].imshow(np.abs(kspace_data), cmap='gray', vmax=0.005*np.abs(kspace_data).max())
        
        # 显示重建的图像
        axs[i, 1].imshow(np.abs(image_rec0), cmap='Greys_r')

        axs[i, j].set_xlabel(f"nb_samples = {nb_samples}")
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"decay = {decay}")
    

axs[0, 0].set_title("k-space noisy data")
axs[0, 1].set_title("Zero-order recon")

# Add a main title (suptitle) to the figure
fig.suptitle(f"k-space noisy data and its Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask = "VDS_along_parallel_lines", signoise = 10, nb_samples = 64, we plot 7 reconstructed images with decay varying in vector_decay = [0.1, 0.5, 0.8, 1, 1.5, 4., 7.] :

• Initially, as decay increases, the quality of the reconstructed image improves until decay ≃ 0.5, after which the image quality starts to degrade.

• To quantitatively evaluate the reconstruction accuracy, we can calculate the Structural Similarity Index between the true image and the reconstructed image.

$\leadsto$   We can compute the Structural Similarity Index to evaluate the quality of reconstruction :

$$ SSIM(x,y) = \frac{(2m_x m_y + C_1)(2\Sigma_{xy} + C_2)}{(m_x^2 + m_y^2 + C_1)(\Sigma_{xx}^2 + \Sigma_{yy}^2 + C_2)} $$

where $m$ represents the mean, and $\Sigma$ denotes the covariance, and $C_1$ and $C_2$ are constants.

The result represents the overall similarity between the two images, ranging from [0, 1], where 1 indicates they are identical, and 0 means there is no similarity at all.

In Python, import the module :

from skimage.metrics import structural_similarity as ssim

Then, simply use the following line to compute the SSIM score and map :

ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)

from skimage.metrics import structural_similarity as ssim

mask="VDS_along_parallel_lines"
vector_signoise = [0, 10, 20, 100]
vector_decay = np.linspace(0, 10, 200)
vector_nb_samples = [int(img_size / i) for i in [2, 4, 8, 16, 32]]

vector_ssim_score = []
vector_image_rec0 = []

'''
mask == "VDS_along_parallel_lines"
signoise = 10
decay = 1.
nb_samples = (int)(img_size/4): 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
'''

for i in range(len(vector_decay)) :
    

        signoise = vector_signoise[1]
        decay = vector_decay[i]
        nb_samples = vector_nb_samples[1]

        kspace_data, image_rec0 = Reconstruction (mask, signoise, decay, nb_samples)

        vector_image_rec0.append(image_rec0)
        ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)
        vector_ssim_score.append(ssim_score)
    

plt.figure(figsize=(10, 6))
plt.plot(vector_decay, vector_ssim_score, linestyle='-', color='red', marker='.', markersize=8, markerfacecolor='red', label='SSIM Score')

plt.xlabel('Decay')
plt.ylabel('SSIM Score')
plt.title('SSIM Score vs Decay')
plt.grid(True)
plt.legend()
plt.show()

$\leadsto$   Remember, cdf_pdecay is obtained by computing the probability density function (PDF) p_decay = np.power(np.abs(kspace_lines),-decay) , followed by normalization as p_decay = p_decay/np.sum(p_decay) , integrating this PDF cdf_pdecay = np.cumsum(p_decay)

The results show that the maximum SSIM Score occurs at decay = 0.8 . Below, we display the reconstructed images for it :

• Decay = 0.8 (corresponding to the maximum SSIM Score)

Where we added a Gaussian complex-valued random noise with signoise = 10 and fixed the nb_samples = (int)(img_size/4) = 64 .

# Find the index and value of the ；maximum Frobenius norm.
max_index = np.argmax(vector_ssim_score)
image_corresponding_max_index = vector_image_rec0[max_index]
decay_corresponding_max_index = vector_decay[max_index]
# print(decay_corresponding_max_index)

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

axes[0].imshow(np.abs(image_corresponding_max_index), cmap='Greys_r')
axes[0].set_title(f'Decay (for Max SSIM Score) : {decay_corresponding_max_index:.2f}')
axes[0].axis('off')

axes[1].imshow(mri_img, cmap='Greys_r')
axes[1].set_title("True image")
axes[1].axis('off')

plt.tight_layout()
plt.show()

• When we select decay = 0.80 , the best reconstructed image is obtained in terms of the SSIM Score.

Q2  

• Second block: Vary parameter nb_samples to define the number of phase encoding lines and study its impact on image quality.

mask="VDS_along_parallel_lines"
vector_signoise = [0, 10, 20, 100]
vector_decay = [0.1, 0.5, 0.8, 1, 1.5, 4., 7.]
vector_nb_samples = [int(img_size / i) for i in [1, 2, 4, 8, 16, 32]]

'''
mask == "VDS_along_parallel_lines"
signoise = 10
decay = 1.
nb_samples = (int)(img_size/4): 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
'''
fig, axs = plt.subplots(len(vector_nb_samples), 2, figsize=(24, 24))

for i in range(len(vector_nb_samples)) :
    for j in range(2):

        signoise = vector_signoise[1]
        decay = vector_decay[3]
        nb_samples = vector_nb_samples[i]

        kspace_data, image_rec0 = Reconstruction (mask, signoise, decay, nb_samples)

        # 显示被噪声影响的 k-space 数据
        axs[i, 0].imshow(np.abs(kspace_data), cmap='gray', vmax=0.005*np.abs(kspace_data).max())
        
        # 显示重建的图像
        axs[i, 1].imshow(np.abs(image_rec0), cmap='Greys_r')

        axs[i, j].set_xlabel(f"nb_samples = {nb_samples}")
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"decay = {decay}")
    

axs[0, 0].set_title("k-space noisy data")
axs[0, 1].set_title("Zero-order recon")

# Add a main title (suptitle) to the figure
fig.suptitle(f"k-space noisy data and its Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask = "VDS_along_parallel_lines", signoise = 10, decay = 1. , we plot 6 reconstructed images with nb_samples varying in vector_nb_samples = [int(img_size / i) for i in [1, 2, 4, 8, 16, 32]] :

• As nb_samples decreases, the quality of the reconstructed image degardes.

• To quantitatively evaluate the reconstruction accuracy, we can calculate the Structural Similarity Index between the true image and the reconstructed image.

$\leadsto$   We can compute the Structural Similarity Index to evaluate the quality of reconstruction :

$$ SSIM(x,y) = \frac{(2m_x m_y + C_1)(2\Sigma_{xy} + C_2)}{(m_x^2 + m_y^2 + C_1)(\Sigma_{xx}^2 + \Sigma_{yy}^2 + C_2)} $$

where $m$ represents the mean, and $\Sigma$ denotes the covariance, and $C_1$ and $C_2$ are constants.

The result represents the overall similarity between the two images, ranging from [0, 1], where 1 indicates they are identical, and 0 means there is no similarity at all.

In Python, import the module :

from skimage.metrics import structural_similarity as ssim

Then, simply use the following line to compute the SSIM score and map :

ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)

from skimage.metrics import structural_similarity as ssim

mask="VDS_along_parallel_lines"
vector_signoise = [0, 10, 20, 100]
vector_decay = [0.1, 0.5, 0.8, 1, 1.5, 4., 7.]
vector_nb_samples = [int(i) for i in np.linspace(1, 256, 256)]

vector_ssim_score = []
vector_image_rec0 = []

'''
mask == "VDS_along_parallel_lines"
signoise = 10
decay = 1.
nb_samples = (int)(img_size/4): 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
'''

for i in range(len(vector_nb_samples)) :
    

        signoise = vector_signoise[1]
        decay = vector_decay[3]
        nb_samples = vector_nb_samples[i]

        kspace_data, image_rec0 = Reconstruction (mask, signoise, decay, nb_samples)

        vector_image_rec0.append(image_rec0)
        ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)
        vector_ssim_score.append(ssim_score)
    

plt.figure(figsize=(10, 6))
plt.plot(vector_nb_samples, vector_ssim_score, linestyle='-', color='red', marker='.', markersize=8, markerfacecolor='red', label='SSIM Score')

plt.xlabel('NB Samples')
plt.ylabel('SSIM Score')
plt.title('SSIM Score vs NB Samples')
plt.grid(True)
plt.legend()
plt.show()

$\leadsto$   Remember, samples is a set of nb_samples uniformly distributed random numbers generated within the range [0, 1] : samples = ra.uniform(0, 1, nb_samples)

The results show that the SSIM score increases as nb_samples increases. Below, we display the reconstructed image for the maximum SSIM score:

• nb_samples = 251 (corresponding to the maximum SSIM Score)

Where we added a Gaussian complex-valued random noise with signoise = 10 and fixed the decay = 1. .

# Find the index and value of the ；maximum Frobenius norm.
max_index = np.argmax(vector_ssim_score)
image_corresponding_max_index = vector_image_rec0[max_index]
nb_samples_corresponding_max_index = vector_nb_samples[max_index]
# print(nb_samples_corresponding_max_index)

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

axes[0].imshow(np.abs(image_corresponding_max_index), cmap='Greys_r')
axes[0].set_title(f'NB Samples (for Max SSIM Score) : {nb_samples_corresponding_max_index:.2f}')
axes[0].axis('off')

axes[1].imshow(mri_img, cmap='Greys_r')
axes[1].set_title("True image")
axes[1].axis('off')

plt.tight_layout()
plt.show()

• When we select nb_samples = 251 , the best reconstructed image is obtained in terms of the SSIM Score.

Code given by Github

import numpy as np
import sys
import numpy.random as ra

eps = sys.float_info.epsilon     # sys.float_info.epsilon 是 Python 中 浮点数最小可区分值（大约为 2.22e-16）。

# 0. 初始化一个 k-space 采样掩码矩阵
'''
c 是一个 形状为 (1, img_size) 的 全 1 的数组（数据类型为 float64）。
np.tile 作用是将数组 c 沿垂直方向（行）复制 img_size 次，生成一个 大小为 (img_size, img_size) 的全 1 矩阵。

'''
c = np.ones((1, img_size), dtype="float64")     
kspace_mask = np.tile(c, (img_size, 1)) 

# 1. 生成伪随机的采样密度和 CDF（累计分布函数）
#kspace_lines = np.linspace(-1/2., 1/2.,img_size)
'''
kspace_lines：生成从 -img_size/2 到 img_size/2 的线性空间，表示 k-space 的相位编码方向上的采样线。
p_decay: 定义一个采样密度分布，它随距离中心频率的距离增大而衰减（decay 控制衰减速度）。
cdf_pdecay: 计算采样密度的累计分布函数 (CDF)，用于将均匀随机变量映射到指定的采样分布中。
'''
kspace_lines = np.linspace(-1/2., 1/2.,img_size)*img_size   
# define the taret sampling density (ie non-iniform over k-space lines)
decay = 1.
# Define the sampling density
p_decay = np.power(np.abs(kspace_lines),-decay)     
p_decay = p_decay/np.sum(p_decay)
# generate its CDF
cdf_pdecay = np.cumsum(p_decay)

pmax = p_decay.max()
pmin = p_decay.min()

# Plot the density and its cumulative distribution function (CDF) p_decay
plt.figure()
fig, axs = plt.subplots(1, 2, figsize= (5,5) )
axs[0].plot(kspace_lines, p_decay)
axs[1].plot(kspace_lines, cdf_pdecay)

# 2. 生成伪随机采样线 
# Perform pseudo-random sampling: technique used:
# draw uniform variables and invert the CDF to get back to p_decay-distributed sampled
nb_samples = (int)(img_size/4)              # 选择采样线的总数，约为总 k-space 线数的 1/4。 256/4 = 64
print(nb_samples)
samples = ra.uniform(0, 1, nb_samples)      # 在 [0, 1] 范围内生成 nb_samples 个均匀分布的随机数。 
                                            # 每个随机数的值都在 0 到 1 之间，且每个随机数的概率是相等的。
'''
np.argwhere(cdf_pdecay == min(cdf_pdecay[(cdf_pdecay - r) > 0]))：
通过反转 CDF，将均匀随机变量 samples 映射到具有指定分布的采样线 gen_klines。
'''
'''
for r in samples: 遍历随机样本 r（r 是 0 到 1 之间的随机数）
cdf_pdecay[(cdf_pdecay - r) > 0]: 找到 CDF 中第一个大于随机数 r 的位置
np.argwhere(cdf_pdecay == min(cdf_pdecay[(cdf_pdecay - r) > 0])): np.argwhere 找到这个最小值在 cdf_pdecay 中的位置（索引）
kspace_lines[np.argwhere(...)]: 从 kspace_lines 中取出对应的采样线位置
gen_klines = [int ]: 将结果转换为整数并存入列表 gen_klines

-> 将随机样本 samples 映射到 cdf_pdecay，并从 kspace_lines 中选择相应的采样线，生成一个伪随机采样线列表 gen_klines。
-> 这些采样线遵循 p_decay 定义的密度分布，低频采样较多，高频采样较少。
'''
gen_klines = [int(kspace_lines[np.argwhere(cdf_pdecay == min(cdf_pdecay[(cdf_pdecay - r) > 0]))]) for r in samples]

# 3. 对采样线进行偏移和平移
# shift the samples lines by half of the k-space (ie image size as we're in Cartesian ref)
'''
生成的采样线 gen_klines 在 [-img_size/2, img_size/2] 范围内，
需要平移到 [0, img_size] 以匹配图像矩阵的索引格式。
'''
gen_klinesb = ((np.array(gen_klines) - 1) / 1).astype(int) + (int)(img_size/2)
#gen_klinesb = ((np.array(gen_klines) - 1) / 1).astype(int)
print(gen_klinesb)   # unsorted samples

# 4. 绘制采样密度分布和直方图
times = np.arange(1, img_size, 1)
# np.bincount(gen_klinesb) 计算每条采样线的出现次数，得到伪随机采样线的实际分布。
lc = np.bincount(gen_klinesb, minlength=len(times))     

# 绘制采样线分布直方图并与原始密度分布 p_decay 比较，检查采样分布是否与目标分布一致。
# check that histogram of sample values fits the prescribed density p_decay
plt.figure()
plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Assigned k-space lines')
plot2, = plt.plot(p_decay,'g',label='Original Spectrum')
plt.xlabel('k-lines')
plt.ylabel('Probability')
plt.legend(handles=[plot1,plot2])
plt.show()

#print(p_decay.min())
sampled_klines = np.array(np.unique(gen_klinesb))       # np.unique() 用于去除重复值，返回唯一的采样线索引列表。
print(sampled_klines)

nblines = np.size(sampled_klines)
print(nblines)
threshold = 2. * p_decay.min()  # sys.float_info.epsilon \simeq 2e-16
kspace_mask = np.zeros((img_size,img_size), dtype="float64")
kspace_mask[sampled_klines,:] = np.ones((nblines,img_size) , dtype="float64")

plt.figure()
plt.imshow(kspace_mask, cmap='gray')
plt.show()

64
[127 127  59 124 137 127 127 127 121 102 127 121 183 126 126 132 124 225
  16 159 128 130 132 123  73 118  78 230 129 116 122 127  85 127 128 127
 127 127 118 127 147 121 127 142 208 126 127 163  30 146 127 127 105  88
 130 127 151 128 127 127 127 114 126 133]

<Figure size 640x480 with 0 Axes>

[ 16  30  59  73  78  85  88 102 105 114 116 118 121 122 123 124 126 127
 128 129 130 132 133 137 142 146 147 151 159 163 183 208 225 230]
34

#import numpy.fft as fft
norm = "ortho"
#norm = None
def fft(x):
    return np.fft.fft2(x, norm=norm)

def ifft(x):
    return np.fft.ifft2(x, norm=norm)


# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise

# Mask data to perform subsampling
kspace_data *= kspace_mask

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))

fig, axs = plt.subplots(2, 2, figsize=(10, 10) )
axs[0,0].imshow(mri_img, cmap='Greys_r')
axs[0,0].set_title("True image")
axs[0,1].imshow(kspace_mask, cmap='gray')
axs[0,1].set_title("Sampling mask")
axs[1,0].imshow(np.abs(kspace_data),  cmap='gray', vmax=0.005*np.abs(kspace_data).max())
#axs[1].imshow(np.abs(np.fft.ifftshift(kspace_data)), cmap='Greys_r')
axs[1,0].set_title("k-space noisy data (1D VDS)")
axs[1,1].imshow(np.abs(image_rec0), cmap='Greys_r')
axs[1,1].set_title("Zero-order recon")
plt.show()